# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

from fractions import Fraction
from numbers import Number
import math

def in_ideal(f, generators):
    """Determine whether a polynomial is in an ideal, or not."""
    g = find_principal(generators)
    if zero_polynomial(divide_polynomials(f, g)[1]):
        return True
    else:
        return False

def find_principal(generators):
    """Find principal ideal, given a list of generating polynomials.

    Arguments can either be a single polynomial, a single element list
    containing a polynomial, or a list of polynomials."""
    if is_polynomial(generators):
        return generators
    elif len(generators) == 1:
        return generators[0]
    else:
        return gcd(generators[0], generators[1:])

def gcd(f, rest):
    """Find the greatest common divisor of a list of polynomials."""
    # Calculuates the gcd recursively.
    if len(rest) > 1 and isinstance(rest[0], list):
        return gcd(f, gcd(rest[0], rest[1:]))
    else:
        # The following two lines of code pulls a single polynomial out of a 
        # single element list.
        # Eg. [[Fraction(-1), Fraction(2)]] => [Fraction(-1), Fraction(2)]
        if isinstance(rest, list) and isinstance(rest[0], list):
            rest = rest[0]
        h = f
        s = rest
        while not zero_polynomial(s):
            remainder = divide_polynomials(h, s)[1]
            h = s
            s = remainder

    return h            

def divide_polynomials(f, g):
    """Divide two polynomials. Returns a tuple of the result and remainder """
    if zero_polynomial(g):
        return (None, f)

    divisor = [0]*max(degree(f), degree(g))
    rem = f

    while not zero_polynomial(rem) and degree(g) <= degree(rem):
        divisor = add_polynomials(
            divisor, monomial_to_polynomial(
                divide_monomials(leading_term(rem), leading_term(g))))

        rem = add_polynomials(
            rem, map(lambda x: -x, multiply_polynomials(
                    monomial_to_polynomial(
                        divide_monomials(
                            leading_term(rem), leading_term(g))), g)))

    return (divisor, rem)


def leading_term(f):
    """Find the leading term of a polynomial."""
    for i in range(len(f)):
        if f[len(f)-i-1] != 0:
            return f[len(f)-i-1], len(f)-i-1

    return 0, float("-inf")

def degree(f):
    """Find the degree of a polynomial."""
    return leading_term(f)[1]

def add_polynomials(f, g):
    """Add two polynomials."""
    h = [0]*max(len(f), len(g))

    for p in [f, g]:
        for i in range(len(p)):
            h[i] += p[i]

    h = reduce_polynomial(h)

    return h

def multiply_polynomials(f, g):
    """Multiply two polynomials."""
    fg = [0]*(len(f)+len(g)-1)

    for i_1,v_1 in enumerate(f):
        for i_2,v_2 in enumerate(g):
            fg[i_1+i_2] += v_1*v_2

    return fg

def reduce_polynomial(f):
    """Remove leading zeros, resulting from subtraction."""
    trim = leading_term(f)[1]

    if math.isinf(trim):
        return [0]

    return f[0:trim]

def zero_polynomial(f):
    """Check if all coefficients are zero."""
    for x in f:
        if x != 0:
            return False

    return True

def monomial_to_polynomial(m):
    """Return a polynomial, given tuple containing coefficient and degree"""
    p = [0]*(m[1]+1)
    p[m[1]] = m[0]

    return p

def divide_monomials(m, n):
    """Divide two monomials."""
    deg = m[1] - n[1]
    c = m[0] / n[0]

    return c, deg

def is_polynomial(f):
    """Test whether a given list is a polynomial, or not."""
    for x in f:
        if not isinstance(x, Number):
            return False

    return True

if __name__ == "__main__":
    # Using example from pg. 42 of Ideals, Varieties, and Algorithms
    # x^6 - 1
    f = [Fraction(-1), Fraction(0), Fraction(0), Fraction(0), Fraction(0), Fraction(0), Fraction(1)]
    # x^4 - 1
    g = [Fraction(-1), Fraction(0), Fraction(0), Fraction(0), Fraction(1)]
    # x^2 - 1
    h = [Fraction(-1), Fraction(0), Fraction(1)]

    print in_ideal(h, [f, g])
